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' ABSTRACT 

oo _ 

, In this paper, we extend our previous analysis (|Lodato & R,icell2004l) of the trans- 

port properties induced by gravitational instabilities in cooling, gaseous accretion discs 
to the case where the disc mass is comparable to the central object. In order to do so, 
we have performed global, three-dimensional smoothed particle hydrodynamics simu- 
lations of massive discs. These new simulations show a much more complex temporal 
evolution with respect to the less massive case. Whereas in the low disc mass case 
a self-regulated, marginally stable state (characterized by an approximately constant 
' radial profile of the stability parameter Q) is easily established, in the high disc mass 

l/^ , case we observe the development of an initial transient and subsequent settling down 

' in a self-regulated state in some simulations, or a series or recurrent spiral episodes, 

I with low azimuthal wave number m, in others. Accretion in this last case can therefore 

Q^. be a highly variable process. On the other hand, we find that the secular evolution 

I ' of the disc is relatively slow. In fact, the time-average of the stress induced by self- 

O I gravity results in accretion time-scales much longer than the dynamical timescale, in 

contrast with previous isothermal simulations of massive accretion discs. We have also 
^ ' compared the resulting stress tensor with the expectations based on a local theory of 

transport, finding no significant evidence for global wave energy transport. 
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1 INTRODUCTION 

It is now becoming progressively more clear (both from the 
theoretical and from the observational point of view) that 
gravitational instabilities might play an important role in 
the dynamics of accretion discs, both on the large scales of 
Active Galactic Nuclei (AGN) and on the small scales of 
protostellar and protoplanetary discs. 

At the AGN scale, a powerful observational tool is 
provided by water maser observations of the gas mo- 
tions at distances of ~ 1 pc from the central black hole 
jGreenhiU fc GwimJiQQilKondratko et al.l2005h . which of- 
ten show that the rotation curve can depart significantly 
from a Keplerian profile, suggesting that the disc itself might 
give a significant contribution to the gravitational field. In- 
deed, in the case of the Seyfert galaxy NGC 1068, simple self- 
gravitating disc models are a ble to reproduce the obs erved 
non-Keplerian rotation curve llLodato fc BertinlbOOSah . 

The case for self-gravitating accretion discs is maybe 
even stronger in the c ontext of star and pla net formation. It 
has long been argued jLin fc Pringlelll987l) that in the early 



stages of star formation the disc mass can be sufficiently 
high for the effect o f its self - gravit y to become important. 
Simple estimates by iLarsorJ il984) lead to the conclusion 
that a gravitationally unstable disc is a likely outcome of the 
isothermal collapse of protostellar clouds. The recent discov- 
ery of a 100 Mq disc s urrounding a mass ive (20 M0) proto- 
star in the M17 nebula llChini et alJ2004l would clearly indi- 
cate, if the mass estimates are confirmed, that the formation 
of a massive star undergoes a phase of self-gravitating disc 
accretion. In the case of low mass (proto)stars, disc mass es- 
timates are uncertain, but discs in the upper end of the mass 
distribution for T Tauri stars iLaunhardt fc Sargentll200ll : 
iNatta et all20o3) can have a mass ~ 30% of the central star. 
In addition, there have also been theoretical suggestions that 
the disc self-gravity might be important in the outer disc of 
FU Orionis objects lEeU fc LirJIlOoi : lArmitaee et alJl200ll : 
iLodato fc BertiDll2003bi) ! 

Recently, the disc self-gravity has been considered es- 
pecially in view of the possibility of forming smaller mass 
companions by the direct fragmentation of a gravitationally 
unstable disc, with obvious implications for the formation 
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of g as Riant planets in a protos tellar disc (see, for exam- 
ple, iBosj 120021: iRice et alJl2003l or, for simil ar arguments 
on the galactic scale, looo^nan fc Tanll2004 ). However, it 
has now become clear that the outcome of gravitational in- 
stabilities strongly depends on the thermodynamics of the 
disc. In particular, the fragmentation of a gravitationally 
unstabl e disc rec^uires^ that the disc be able to cool very effi- 
ciently llGamm ic 2001': ' Rice et all2003all . with a cooling rate 
tcooi ^ 3f2~^ (where Q, is the angular velocity of the disc). 
This is a rather severe requirement for a protostellar disc, so 
that the formation of planets from the direct fragmentation 
of a massive disc seems rather unlikely (however, the inclu- 
sion of the effects of conve ction might substantially reduce 
the cooling timescale, see lBosjl2004all . In any case, even 
the simple development of a gravitational instability in the 
form of a spiral structure can significantly accelerate the pro- 
cess of planet formation by favouring the agglomeration o f 
planetesimals llHaghighipour fc Bossi200jl : lRice et alj2004l) . 
Moreover, gravitational stresses might play another impor- 
tant role in the process of planet formation, since the mixing 
induced by these stresses can influence the early evolution 
of chemical elements in protoplanetary discs ^Boss 2004b). 

Another important issue related to the disc self-gravity 
is its ability to transport angular momentum and thus to 
promote the process of accretion. This has been recog- 
nized in the analysis of isothermal simulatio ns of massive 
discs iLaughlin fc Bodcnheimc r 1994: Laughli n fc Rozvczkal 
ll99(J) and subsequently analysed in a num ber of papers 
JPickett et alll998l.l200ol:lNelson et al.l2000h . However, only 
recently have simulations been performed which include a 
relatively more complex tr e atment of the thermal status of 
the d isc llRice et al.ll2003al : iPickett et alj|2003t iMeiia et alJ 
I2OO4I) . In a previous paper llLodato fc Ricdl2004l hereafter 
Paper I) , we have explored in detail the transport properties 
of self-gravitating, cooling discs, by performing some high- 
resolution smoothed particle hydrodynamics (SPH) simula- 
tions. These simulations were used to explore the extent to 
which angular momentum and energy transport via gravita- 
tional instability can be regarded as a local phenomenon, or 
whether the intrinsic global nature of the instabilities would 
preclude such an approach. The simulations described in 
iPaper J included disc masses up to 0.25M* (where is 
the mass of the central object) and showed how, when the 
disc is allowed to cool down due to external cooling and to 
heat up as a consequence of the development of gravitational 
spiral instabilities, it eventually settles down in a marginally 
stable state, characterized by an approximately flat profile 
of the parameter Q: 



Q 



TvGT, 



(1) 



where Cs is the thermal speed, k is the epicyclic frequency 
and E is the surface density of the disc. The resulting spi- 
ral structure turned out to be a quasi-stationary structure, 
evolving only on the long "viscous" time-scale, and able to 
transport angular momentum and energy at a rate which is 
in reasonable agreement with the expectations based on a 
local (viscous) theory of transport. 

It can be shown that, for self-regulated discs, a sim- 
ple relationship holds between the disc aspect ratio H/R 
and the mass ratio M^isc/Mt. The disc is in fact marginally 
stable when its temperature is small enough that H/R < 



MdiBc/Mf Of course, this relationship can only hold in the 
limit where Mdisc/Af, <C 1. As Mdisc/Af, grows, the disc 
becomes thicker, and global effects are going to be more im- 
portant. This was already shown by us in Paper I. The limit 
Mdiac/M, ~ 1 is therefore particularly interesting, because 
we expect a different physics to come into play. 

In this paper we extend our previous analysis, consid- 
ering the case where the disc mass is particularly high. We 
have in fact considered the two cases where Mdisc = 0.5M* 
and A/disc ~ IM*. The results of these new simulations 
present significant differences with respect to our previous 
analysis. In the high disc mass case, we now find the de- 
velopment of transient, global spiral structures, which last 
for roughly one dynamical timescale and develop recur- 
rently during the course of the simulation. While in the 
Afdisc = 0.5Mt case, after a flrst transient, the disc is even- 
tually able to settle down in a quasi-stationary state, in the 
equal mass case the disc in not able to "self-regulate" the 
strength of the spiral disturbance in order to counterbalance 
exactly the externally imposed cooling and short-lived global 
disturbances recurrently develop during the whole simula- 
tion. Accretion in this last case will therefore appear to be 
a highly variable process. 

The paper is organized as follows. In Section|21we briefly 
review the relevant physical concepts of the problem and the 
details of the numerical setup that we adopt. In Section |3 
we describe the results of our simulations. In Section |1| we 
discuss our results and draw our conclusions. 



2 BASIC CONCEPTS AND NUMERICAL 
SETUP 

2.1 Basic physical quantities 

In this section we briefly review the basic concepts about the 
evolution of viscous discs and the development of gravita- 
tional instabilities and their associated tr ansport p roperties. 
More details can be found in section 2 of iPaoerJ and refer- 
ences therein. 

The equations of motion for an axisymmetric disc, in 
cylindrical coordinates, are the equation of continuity: 

rlT 1 f) 

a^ + ^aS(^^-) = °' (2) 

and the azimuthal component of Euler's equation (expressed 
in terms of angular momentum conservation): 

where u is the radial velocity and Tr^ is the relevant com- 
ponent of the viscous stress tensor, integrated in the vertical 
direction. This last term is the basic ingredient in the the- 
ory of accretion discs. Standard hydrodynamical viscosity 
is not sufficient to provide accretion at the required rates, 
and therefore Tr^ is generally descri bed by means of ad 
hoc p rescriptions. The a-prescription dShakura fc SunvaevI 
Il973r . based on simple arguments on the relevant physical 
scales of turbulent cells in the discs, assumes that Tr^ is 
proportional to the disc pressure: 



dlnn 



dlnR 



(4) 
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where the proportionality constant a is an unknown pa- 
rameter, usually considered to be smaller than unity. The 
a-prescription can be also put in the following equivalent 
form, which involves the kinematical viscosity coefficient u: 



(5) 



where H is the thickness of the disc. 

The effect of viscosity on the energy balance is twofold: 
viscous torques convect energy between neighbouring annuli 
of the disc and they dissipate energy. The energy which is 
convected across a ring at radius R per unit time is given 
by: 
H W 

— = 2nR''TR4,n, (6) 

while the dissipation rate per unit surface is given by: 

D{R)=Tr4,\RQ!\. (7) 

If the disc is in thermal equilibrium, we can derive a useful 
relation between the viscosity coefficient a and the cooling 
timescale. If we assume that cooling can be simply parame- 
terised in the following way: 

Q--^-..rM^, (8) 



7(7 - l)fcool ' 



where U is the internal energy per unit surface and 7 is 
the ratio of the specific heats, then, equating the viscous 
dissipation term, expressed in Eq. Q to the cooling term, 
expressed in Eq. ||HJ, leads to: 



dlni? 



7(7 - l)tcoolfi' 



(9) 



where we have used also Eq. 0. 



2.2 Diagnostics for gravitationally induced 
transport 

In order to describe the development and the properties of 
a self-gravi tating st ructure, we will use the same diagnostics 
adopted in Paper I, and we summarize them here. 

The basic parameter that describes the local, axisym- 
metric stability of a gaseous disc with respect to gravita- 
tional instabilities is the Q parameter defined in Eq. Q. It 
actually turns out that the stability criterion based on this 
parameter, namely that Q > 1, is quite robust, and gener- 
ally applicable also in the case of global, non-axisymmetric 
disturbances (see also Paper I). As mentioned in the Intro- 
duction, and as discussed in more details elsewhere (Paper I, 
^ertin 1997; Paczvriski 1978} the competitive effects of ex- 
ternal cooling and effective heating due to gravitational in- 
stabilities tend to result in a self-regulated state, where the 
value of Q is constant throughout the disc and very close 
to its mar ginal stability value. The simulations presented 
in IPaper ll confirmed the effectiveness of the self-regulation 
mechanism. We will therefore track the development of a 
self-regulated state by looking at the profile of Q. 

The relevant component of the stress tensor as- 
sociated with gr avitational instabilities is defined as 
jLvnden-Bell fc Kalnai&.197a'l : 



Az 
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where gn and are the radial and azimuthal component 
of the gravitational field, respectively. The full stress tensor 
also includes the 'Reynolds' stress, defined as: 



(11) 



where (5v = v — u, with v the fluid velocity and u the mean 
fluid velocity. It is always possible to deflne an effective a 
parameter associated with gravitational instabilities, by set- 
ting: 



r, 
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A local description of angular momentum transport then 
requires that the stress tensor, and therefore a, at a given 
location is only determined by local conditions in the disc. 
In Paper I we have already shown that in self-gravitating 
discs a region of size ~ 10-ff contributes to the stress at a 
given location, therefore implying that for discs as thick as 
H/R ~ 0.1 the stress tensor is not determined locally. 

In addition, there is also another issue concerning the 
use of a viscous formalism in self-gravitating accretion discs. 
In fact,_Balbus & Papaloizou ( 1992i) have shown that in gen- 
eral the energy transport provided by gravitational insta- 
bilities contains global terms, associated with wave energy 
transport, which cannot be directly associated with an ef- 
fective viscosity. It remains to be discussed how large these 
term s are, and under which conditions they play a signiflcant 
role. iBalbus fc Papaloizoul il999l) argue that global energy 
transport should not play an important role in a marginally 
stable disc. In this case, a viscous description of transport 
should be appropriate and therefore energy should be trans- 
ported according to Eq. Q and it should be dissipated ac- 
cording to Eq. (|7J . If the disc is in thermal equilibrium, the 
latter statement is equivalent to the prescription that a and 
tcooi are related through Eq. Q . This offers us a simple test 
to check whether energy dissipation is actually viscous or not 
in massive discs. In particular, since in all our simulations 
(as discussed below), we adopt 7 — 5/3 and tcooi = 7.511"^, 
the expected value of a for viscous dissipation to balance 
cooling is a « 0.053. 

2.3 Numerical issues and setup 

The numerical setup is largely identical to that used in Pa- 
per I. We use Smoothed Part icle Hydrodynamics (SPH) (see 
lBendll990l: lMonaghanlll992r) to perform three-dimensional 
simulations of gaseous accretion discs. We consider a sys- 
tem comprising a central star, with mass M*, surrounded 
by a disc, with mass Mdisc, and we consider mass ratios, 
q — Mdisc /Mt, of 0.5 and 1. The disc is modelled using 
250000 SPH particles while the central object is modelled 
as a point mass onto which gas particles may accrete if they 
move to within the accretion radius. Both the gas particles 
and the point mass use a tree to determine n eighbours and 
to calculate gravitational forces jBendll990D . and the cen- 
tral object is free to move under the gravitational influence 
of the disc gas. 

The disc particles are distributed such as to give an ini- 
tial surface density proflle of E oc R"^ with an initial tem- 
perature profile chosen to be T cx R~^^^. The initial velocity 
profile is calculated by including the enclosed mass and the 
pressure gradient when determining the angular frequency 
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Q.. With these initial conditions, the minimum value of Q is 
attained at the outer edge of the disc. For each simulation, 
the temperature normalisation is chosen such that the min- 
imum value of Q is Qmin = 2, giving a disc that is initially 
gravitationally stable. The disc is also assume d to be in vert i- 
cal hydrostatic equilibrium (see, for example. IPringljll981^ . 
The particles are therefore distributed such that the vertical 
density profile is a Gaussian with a scaleheight H given by 
H = Ca/n. Actually, in a self-gravitating disc, the vertical 
density profile is not rigorously Gaussian jBertin fc Lodatol 
ll999l:lHurdl2OO0ll . so our initial setup, strictly speaking, is 
not in dynamical equilibrium. However, dynamical equilib- 
rium is achieved rapidly (i.e., in a dynamical timescale) dur- 
ing the simulation. 

Our calculations are essentially scale free. In dimension- 
less units, the disc extends from _Rin = 0.25 to iiout = 25 
and the star has a mass of M« = 1. In these units, the orbital 
period at i? = 1 is 27r code units and one orbital period at 
the outer edge of the disc is approximately 800 time units. 

Since the main aim of this work, as in Paper I, is to 
investigate transport properties of self-gravitating accretion 
discs, we explicitly solve the energy balance for the gas. We 
allow the disc to heat up due to both PAV work and viscous 
dissipation and assume that the ratio of the spec ific heats 
is 7 = 5/3. We then impose a cooling of the form iGammid 
I2OOII) 

OL t'cool 

where Ui is the internal energy per unit m ass and tcooi is 
given by tcooi = ■ Previous simulations taammi3l200ll: 

iRice et al. 20 03b'l have shown that discs fragment for /3 ^ 3 
with this fragmentat ion boundary poss ibly increasing as the 
disc mass increases iRice et alj|2003bl) . To ensure that our 
discs do not fragment, we use (3 = 7.5. The imposed cool- 
ing should lead to the growth of the gravitational instability 
which will heat the disc, returning it to a state of marginal 
stability. There will also be additional heating through the 
artificial SPH viscosity. I n these simulation we use the stan- 
dard SPH viscositv fe.g. |Monagha nlll992ll but also include 
the Balsara switch (lBalsaral ^9oi~to reduce the shear com- 
ponent of the artificial viscosity. As discussed in detail in 
Paper I, the angular momentum transport induced by the 
artificial viscosity is at least a factor of 10 smaller than that 
due to the gravitational instabilities. The artificial viscosity 
should therefore play only a minor role in the disc dynamics. 



3 RESULTS 

In this section we describe the main results of our simula- 
tions. As mentioned above, all the simulations are scale-free 
and can therefore describe the basic dynamics of different as- 
trophysical systems (from protostellar discs, to AGN discs 
or even to galactic discs), even if we do not intend to model 
the detailed physical appearance of a specific one. However, 
when numbers are needed, we will consider as a reference 
case, for illustrative purposes, a "low mass proto-star cir- 
cumstellar disc" , where the central object mass is IMq and 
the unit length of the simulations is J?q = 1 au. Thus, in 
this reference case, the disc extends from _Rin = 0.25 au to 
RovLt ~ 25 au. We will express all times in units of the Kep- 



lerian rotation period of the outer disc to ~ 27r/r2out, where 
^^out = GMi,/Rauf For most cases, we have followed the 
simulations up to « 7to, i.e. up to a few thermal timescales 
at the outer edge of the disc. 



3.1 The Mdisc = 0.5A/* case 

The evolution of the disc in this case is characterized by two 
phases: (?) soon after the beginning of the simulation the disc 
develops a large-scale, transient spiral structure that subse- 
quently vanishes, (ii) at later times the disc settles down 
in a quasi-stationary state, characterized by a much slower 
evolution. Here we describe the two stages separately. 

3.1.1 The initial transient 

In the earliest stages of the simulation the disc simply cools 
down, until (at roughly t = 2to) the minimum value of Q ap- 
proaches unity and the first gravitational disturbances start 
to develop. These are dominated by a grand-design, m — 2 
spiral mode. The pattern frequency of this mode (i.e. the 
angular velocity of a reference frame in which the spiral is 
stationary) is Qp « 1.5nout. In this case the corotation ra- 
dius (i.e. the radius at which the spiral mode and the gas 
rotate with the same angular velocity) is i? ~ 19, very close 
to the outer edge of the disc. This large-scale structure is 
however not long-lived, and has essentially disappeared at 
t ~ 3to (i.e. after one dynamical timescale), leaving behind a 
smaller-scale disturbance characterized by modes with larger 
m and confined within the inner disc. Fig.^shows the loga- 
rithm of the disc surface density (with a colour scale ranging 
between 10* — 10^ g cm^^) during the development of this 
transient disturbance. The four images refer tot — Q, 1.9, 2.5 
and 3to, respectively. The linear scale of the images ranges 
between -25 and 25 for both axes. 

The development of this spiral is accompanied by a sig- 
nificant modification of the profile of Q. Fig. |5| shows the 
Q-profile at t = l.Qto (solid line), t = 2.7to (dotted line) 
and, for comparison, at t — 7to (i.e. at the end of the sim- 
ulation, when a self-regulated state is eventually achieved; 
dashed line). It can be seen that this first spiral transient 
leads to a strong "heating" of the outer disc, with Q rising 
well above unity. In this way the outer disc is stabilized and 
the large-scale spiral structure rapidly vanishes. 

Another important feature associated with this tran- 
sient spiral instability is the mass transport induced by it. 
We have evaluated the strength of the torques excited by 
the gravitational instability as discussed in section |5| The 
left panel of Fig. 2] shows the evolution of the parameter a 
at four different times during the transient: t — 2.28to (solid 
line), t — 2.52to (dashed line), t — 2.69to (long-dashed line) 
and t — 2.93to (dot-dashed line). The two dotted lines show 
the expected "equilibrium" value a = 0.053, obtained from 
balancing viscous heating and cooling (as discussed above in 
section l^l, and the maximum torque strength (q = 0.005) 
due to artificial SPH viscosity (see Paper I, Appendix). It 
is apparent that the transient induces a strong angular mo- 
mentum transport in the outer disc, with maximum values of 
a raising above 0.2. At later times, when the transient van- 
ishes, the value of a decreases accordingly, becoming much 
smaller than 0.053 in the outer disc, and close to 0.053 within 
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Figure 1. Time evolution of the logarithm of the surface density S (with a colour scale ranging between 10* - 10^ g cm 2) during 
the development of the transient spiral structure at the beginning of the M^isc = 0.5M* simulation. They refer to t = (upper left), 
t ft! 1.9to (upper right), t ^ 2.5to (lower left) and t ^ 3to (lower right). The linear scale of the images ranges from -25 to 25 (in code 
units) for both axes. 



R = 5. The right panel of Fig. |31 shows the corresponding 
mass accretion rates (if we use the scale units of the reference 
case). The AI induced by the transient is significant in the 
outer disc, becoming as large as « 10~'^ Mq/jt, and then 
decreasing significantly once the transient has disappeared. 

The large mass accretion rate of course induces a signif- 
icant redistribution of matter in the disc. This is shown in 
Fig. where we plot the azimuthally averaged surface den- 
sity of the disc ai t — (solid line) and ai t — 2.93to (dotted 
line). The spiral structure has produced a steepening of the 
surface density profile. 

iBertin fc Lodatol 1^9^ have studied a class of steady- 



state, self-regulated massive accretion disc models. The 
surface density profile of these models asymptotically ap- 
proaches S ~ at large radii, but it shows significant de- 
partures from a simple power-law at small radii, becoming 
progressively steeper. These self-regulated models depend on 
only one dimensional parameter: a scale- length Rs, defined 
as follows: 



= 2GM* 




where Q ~ 1 is the marginally stable value of Q. The dashed 
line in Fig. |3 shows the surface density profile of this ana- 
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Figure 2. Evolution of tlie profile of Q during the transient in 
the Af(jist. = O.SAf* simulation. The three curves refer tot = 1.9to 
(solid line), t = 2.7to (dotted line) and t = 7to (dashed line). 




log(R/ au) 

Figure 3. Azimuthally averages surface density profiles at t = 
(solid line) and at t = 2.93to (dotted line). The dashed line shows 
the surface density profile predicted by the self-regulated disc 
models of iBertin fc Lodatoi il999i) . 



lytical model, with Rg — 200 au. With this choice, the mass 
enclosed within 25 au is 0.5M^, i.e. the same as in our simu- 
lation. This density profile matches well the surface density 
profile obtained after the occurrence of the transient spiral. 
This suggest that this episode of gravitational instability is 
the response of the disc to the initial conditions assumed 
here. The disc's "preferred" state is a self-regulated state. 



in which the surface density profile is slightly steeper than 
the initial E cx R~^. The large-scale spiral transient that we 
observe has the effect of moving just enough matter so as to 
steepen the surface density profile to match the one implied 
by the self-regulation condition. 



3.1.2 Long-term evolution 

After the first transient, the value of Q in the outer disc 
is relatively high (see Fig. |5J and the spiral structure is 
confined to the inner disc. In this way, the main source of 
heating in the outer disc (i.e., the effective heating due to 
gravitational instabilities) is turned off and the disc cools 
down until a new generation of gravitational disturbances 
eventually develops. Fig. |K| shows the surface density of the 
disc at the end of the simulation. The scales are exactly the 
same as in Fig. 

This second episode of gravitational instability is how- 
ever different in nature to the previous one. It is not tran- 
sient, but lasts for several outer dynamical timescales (from 
t ~ 4:to until at least the end of the simulation, ai t — 7to), 
showing no significant evolution throughout. This can be 
seen in Fig.|n] which shows the profile of Q at t = 5.2to (solid 
line), t — 5.8to (dotted line) and t = 6.6to (dashed line) (cf., 
for comparison the evolution shown in Fig.|5J. Fig. |7| shows 
the amplitude of the first 15 Fourier components of the den- 
sity structure at the end of the simulations, computed as 
in Paper I. It can be seen that the disc is dominated by 
low-m modes, with m < 5, even if it does not show a clear 
two-armed structure as during the first transient phase. 

In order to discuss the transport properties associated 
with this quasi-stationary structure, we have computed the 
stress tensor averaged over 640 time units, i.e « O.Sfo- The 
results are shown in Figs. |H| and |^ The bottom panel of 
Fig. |5] shows the average stress tensor, as measured by the 
Q parameter (solid line), and the expected value of a ob- 
tained by imposing that the viscous heating rate exactly 
balances the cooling rate (dotted line). There is substantial 
agreement between the two values. The lower dotted line 
shows, as in the left panel of Fig. 0] the maximum contribu- 
tion to the stress coming from artificial SPH viscosity. The 
upper panel of Fig. |H| shows a comparison between the cool- 
ing rate (dashed line) and D{R), the viscous heating rate 
(solid line), defined in Eq. O. Again, there is a substantial 
agreement between the two suggesting that, even for this 
relatively massive disc, a state of thermal equilibrium has 
been reached and that the energy dissipation process is well 
described by a simple viscous model. 

In Fig.|5]we plot the radial profiles of M (upper panel) 
and of t^/t-K (lower panel), where t^, = R^ /v is the vis- 
cous timescale and t^ = is a measure of the dynamical 
timescale. We see that the mass accretion rate acquires an 
almost constant value, consistent with the quasi-stationary 
character of the spiral structure, with M ~ 2 lO~^M0/yr. 
Note that, even if the value of a is not particularly high, 
the mass accretion rate is conspicuous since the disc is rel- 
atively hot. However, despite the high M, the inspection of 
the bottom panel of Fig. |3]shows that accretion occurs on a 
long timescale, more than three orders of magnitude larger 
than the dynamical timescale. 
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Figure 4. Left: Evolution of the gravitational torque strength (as measured by the dimensionless parameter a, during the transient. 
The four lines refer to t = 2.28to (solid line), t = 2.52to (dashed line), t = 2.69to (long-dashed line) and t = 2.93to (dot-dashed line). 
The two dotted lines show the expected "equilibrium" value a = 0.053, obtained from balancing viscous heating and cooling, and the 
maximum torque strength (o = 0.005) due to artificial SPH viscosity (see Paper I, Appendix). Right: Corresponding mass accretion 
rates, assuming that the central object is IMq and that the unit length is 1 au. 





Figure 5. Surface density in the A/disc = 0.5M* case, at the end 
of the simulation. The scales are the same as in Fig. 
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Figure 6. Profiles of Q at t = 5.2to (solid line), t - 
line) and t = 6.6fo (dashed line). 
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3.2 The Afdisc = IM* case 

In this case, the initial evolution of the disc is similar to 
that observed in the other simulations. The disc initially 
simply cools down until Q becomes of the order of unity in 
the inner disc (where the cooling time is fastest) and grav- 
itational instabilities set in. However, the long-term evolu- 



tion of this massive disc differs from the long-term evolu- 
tion of the lower mass discs. A relatively small-amplitude 
spiral structure, confined within R « 15, is always present 
once the disc has cooled down sufficiently. However, a series 
of recurrent, prominent large-scale spiral disturbances occur 
repeatedly throughout the course of the simulation. The am- 
plitude of the m = 2 Fourier component, in the outer disc, 
is plotted as a function of time in Fig. 1101 Major episodes of 
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Figure 7. Amplitude of the first Fourier components of the den- 
sity structure in tlie M^i^^ = 0.5M» case for different radial ranges 
in the disc. 
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Figure 9. Bottom: ratio of the accretion timescale = u/B? 
over the dynamical timescale tx = H"^, for the MiXisc = O.SM* 
case. Top: time averaged mass accretion rate at the end of the 
simulation. 
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Figure 8. Bottom: (solid line) Effective a produced by gravita- 
tional instabilities in the Mdisc = O.SAf* case; (dotted line) value 
of a expected from balancing viscous heating and external cool- 
ing. The lower dotted line, shows the maximum contribution to 
the stress given by artificial viscosity in SPH, as in Fig. |^ Top: 
(solid line) Viscous heating power D{R) compared to (dashed 
line) externally imposed cooling. 



spiral activity can be easily recognized at t « l.Sto, i ~ 4to 
and t ~ 6to. A similar cyclic behaviour, with a modula- 
tion of the amplitude of the dominant m = 2 mode, has 
been found in previous simulations of self-RravitatinR discs 



JSellwood & Carlbere 


ll984HLauehlin & BodenheimeJll994l: 


iLauehlin et al.ll99gD. 


Lauehlin et al.l lll99a) have attributed 



this behaviour to the effect of non-linear coupling between 
different modes. 



3.2.1 The large-scale spiral 

Since all the large-scale spiral disturbances share a similar 
physical appearance, we will describe, as an example, one of 
such episodes, which developed at t ~ 5.8io and has vanished 
by t ~ 6.5to. The structure of the disc at t « 5.9to is shown 
in the left panel of Fig. 1111 (for comparison, the right panel 
shows the structure at t ~ 3to, during a period of low spiral 
activity, see Fig. IIUH . The disc is clearly dominated by a 
large-scale m = 2 mode, with a pattern frequency Sip « 
3f2out, resulting in corotation at 7? ~ 14. Fig. 1121 shows the 
evolution of the profile of Q (left) and of the azimuthally 
averaged surface density E (right) during the development 
of the spiral structure. The curves refer to t ~ 5.4fo (before 
the onset of the instability, solid line), t ~ 5.9to (during 
the disturbance, dotted line), and at t ~ 6.5to (after the 
spiral episode has vanished, dashed line). Before the spiral 
develops, the disc is characterized by a flat, self-regulated, 
Q profile out to i? « 17. Therefore, it appears that the 
corotation of the spiral mode occurs close to the outer edge 
of the self-regulated portion of the disc. Note that the inner 
Lindblad resonance (ILR) for this mode occurs at J? « 5, 
and does not influence the propagation of the density wave 
to small radii, in our hydrodynamical simulations. On the 
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Figure 10. The Mjiat. = M« simulation. Amplitude of the m = 
2 Fourier component of the disc surface density as a function 
of time. Major episodes of spiral activity can be recognized at 
t 1.5*0, t ~ 4to and t « 6to. 

other hand, it is well known that, for a stellar (coUisionless) 
disc, waves are efficiently absorbed at the ILR. 

The effect of the development of the spiral is to remove 
angular momentum from inside corotation and to release it 
to the matter located outside corotation. As a result, outside 
corotation the disc expands and its surface density decreases. 
This is evidenced from the appearance of a local minimum 
in the surface density at corotation (dotted line in the right 
panel of Fig. 1121 , which corresponds to a local maximum in 
the Q profile. After the spiral has vanished, the Q profile is 
very similar to before the onset of the transient structure, 
with the self-regulated part of the disc extending out to a 
slightly smaller radius. 

3.2.2 Secular evolution 

As mentioned above, the appearance of a grand-design m — 
2 mode gives rise to a substantial angular momentum trans- 
port over a small period of time (namely, the duration of the 
m = 2 disturbance). We have computed the time average of 
the the effective a produced by the gravitational instabilities 
at two different times during the evolution of the disc: (i) be- 
tween t — 2.7to and 3.5to, during which period no dominant 
large scale structure was observed, and (ii) between t = 5.8to 
and 6.6to, i.e. during the development of a large scale spiral 
structure (see also Fig. IIUII . The results are shown in Fig. 
1131 where a{R) is shown during low spiral activity (solid 
line) and high spiral activity (dashed line). During high spi- 
ral activity the gravitational torques are larger, but not by 
a particularly significant amount. During both periods the 
value of a lies very close to the expectation based on a bal- 
ance between viscous heating and cooling. If the reference 
numerical values for the relevant scales are used [IMq for 
the unit mass and 1 au for the unit length), the accretion 
rate would be of the order of M ^ 10~*MQ/yr. The accre- 
tion timescale is still very large compared to the dynamical 
time {t^/t-K ~ 10'^), and has a very weak dependence on 
radius, as in the case where Afdisc = 0.5M^. 

It is interesting to compare the torque computed from 
Eq. I12L as shown in Fig. 1131 with the expectations based on 



the linear modal theory of spiral structure in galaxy discs. 
The angular momentum flux associated with a given spiral 
mode can be computed based on the relevant wave action 
and group velocity (see Berlin 198^ . We can thus express 
the angular momentum flux in terms of an effective am, 
associated with a mode with azimuthal wave number m. 
For a tightly wound mode in a light disc we have: 



where A = (5Em/S is the amplitude of the mode, eo = 
-kGTj/Rk^ and k = 2eoRk is the dimensionless radial 
wavenumber of the mode. If we measure the amplitude A 
for the m — 2 spiral mode shown in the left panel of Fig. 
111! and consider the short trailing branch of the dispersion 
relation, outside corotation we obtain am=2 ~ 0.03, which 
shows that a significant fraction of the angular momentum 
transport is indeed associated with the global m = 2 mode. 
Of course some effects beyond the simple linear theory for 
tightly wound modes in light discs, used in the above esti- 
mate, are expected (and should be further analyzed) in the 
present case, in which the structure has reached non-linear 
amplitudes and develops in a relatively heavy disc. 

An important question to be addressed is whether en- 
ergy dissipation during the high spiral activity is intrinsi- 
cally different to standard viscous dissipation, or, in other 
words, whether wave transport of energy is significant in our 
simulations. As discussed in section |5| the test that we have 
generally used in order to assess this issue is a comparison 
between the value of a as computed from the gravitational 
torque strength and the expected value obtained from a bal- 
ance between viscous heating and cooling. In all simulations 
up to now we have interpreted the agreement between the 
two values as a demonstration that wave transport does not 
play a major role. In the present case, with Mdisc = Af*, 
gravitational instabilities are clearly global and different in 
nature with respect to the lower mass cases: here the pattern 
of the instability rapidly changes with time (in particular, 
the strength if the dominant m = 2 mode, see Fig. 1101 and 
the value of Q, although always relatively close to unity, 
never really settles down in a self-regulated state. From Fig. 
I13l we see that even during the high spiral activity, the value 
of a is relatively close to the expected value (dashed line in 
Fig. 1131 . being only slightly larger in the outer disc. It would 
be, however, wrong to interpret this small disagreement as 
an indication that wave energy transport is dominant, since 
the agreement is only expected if the disc is in thermal equi- 
librium. The change in internal energy of the disc as a result 
of the heating due to the development of the spiral can easily 
account for the discrepancy shown in Fig. 1131 

These results do, however, suggest that massive discs 
never settle into a state of thermal equilibrium. If wo in- 
terpret this in the same way as (^aughlin ct al. .199^ . 
non-linear coupling between the different modes results in 
episodes when m = 2 modes heat the outer disc and make it 
more stable, followed be episodes when cooling dominates, 
and the disc is driven back towards a period of instability. 




R 



Figure 12. Evolution of Q (left) and of the azimuthally averaged surface density E (right) during the development of a spiral disturbance 
in the M^^^^ = Ah case. The curves refer to t 5Ato (before the onset of the instability, solid line), t 5.9to (during the disturbance, 
dotted line), and at t 6.34o (after the spiral episode has vanished, dashed line). 



4 CONCLUSIONS 

In this paper we have extended our previous analysis 
of the transport properties of self-gravitating discs, pre- 
sented in Paper I. In particular, here we describe the 
results of two new simulations, in which the disc mass 
is of the same order of the central object mass, namely 
Afdisc = 0.5Af* and A/disc = IM*. Given the specific re- 



gion of the parameter space that we explore, our simula- 
tio ns have two previous analog ues: the A^'-body simulations 
bv lSellwood fc Carlberd il984) . who mimicked the cooling 
of the stellar component of a galaxy disc through continuous 
accret ion of stars, and a series of simulations by Laughlin 
et al. jLauehlin fc Bodenheimeilll994l : iLauehlin et alJll99l 
Il998h . who simulated massive gaseous disc, adopting sim- 
ple equations of state (mostly isothermal) without any spe- 
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cific cooling. The main difference between the latter simu- 
lations and those presented here lies in the way we handle 
the disc thermodynamics. In fact, we explicitly solve the 
energy equation of the disc, and introduce an external cool- 
ing term, simply parameterized in terms of a cooling time 
^cooi = By suitably choosing the parameter f3, we are 

able to prevent the disc from fragmenting, and to produce 
a persistent spiral structure. 

The specific aim of this work is to study the transport 
of angular momentum induced by gravitational instabilities 
and the associated energy dissipation process, to test in par- 
ticular whether wave transp ort of energy plays an impor tant 
role in the disc dynamics jBalbus fc Papaloizoul ITQQgl) . In 
this respect, the results of this paper can be summarized as 
follows; 

(i) The development of a gravitational spiral structure 
is indeed able to transport efficiently angular momentum 
in the disc, hence favouring accretion. The time-averaged 
strength of the the stress induce d by the disc self-gravity, a s 
measured by the parameter a dShakura fc Sunvaevlll973) . 
is a ~ 0.05, so that the associated energy dissipation (com- 
puted from Eq. O) almost balances the externally imposed 
cooling. This shows that global wave transport does not play 
a major role in our simulatio ns, in line with the predictions of 
iBalbus Sz Papaloizod (jlQOO') , who argued that self-regulated 
discs, with Q « 1 (as indeed are ours) should be free from 
wave transport effects. 

(ii) Contrary to the results of iLauehlin fc Bodenheimerl 
il99'J) . the accretion process induced by the disc self-gravity 
occurs on the long 'viscous' time-scale, rather than on the 
much shorter dynamical time-scale, so that a massive self- 
gravitating disc can survive for a relatively long period (see 
Fig-EJ- This important difference can be understood by sim- 
ply considering the parameter Q, as defined in Eq. If 
Q is smaller than unity, gravitational instabilities develop 
and the disc responds so as to increase the value of Q to 
re ach a stable configuration. If t he disc is isothermal (as 
in lLaughlin fc Bodenheimeij *1994'). the only way this can be 
achieved is by significantly reducing E, i.e. by inducing a fast 
accretion process. On the other hand, if the disc is allowed 
to dissipate energy and to heat up, increasing the value of 
Cs (as in our simulations), the disc is stabilized more easily 
and the required amplitude of the gravitational disturbance 
is smaller, hence producing a slower accretion. 

(iii) When the disc mass is small with respect to the cen- 
tral object, as in the case discussed in Paper I, the disc 
rapidly settles down is a self-regulated state where the am- 
plitude of the spiral structure and the radial profile of Q 
change very little with time. In these cases, the spiral struc- 
ture is tightly wound and dominated by high-m modes. In 
the high disc mass case presented here, the evolution is more 
complex. In the case where Mdisc = 0.5AI-^ we observe the 
development of an initial transient m = 2 spiral structure, 
which induces a strong redistribution of matter and a steep- 
ening of the surface density profile, before the disc is eventu- 
ally able to settle down as in the low disc mass case. Similar 
initial transients are also observed byjMeiia et al. ( 2 0o3) . We 
attribute it to an effect of the initial condition adopted. In 
particular, the chosen initial surface density profile E oc 

is actually shallower than analytical esti mates of the sur- 
face density of self-regulated massive discs (iBertin fc Lodatol 
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Figure 13. Time average of the effective a produced by gravi- 
tational instabilities between t = 2.7to and t = 3.5to (where no 
large scale structure is observed; solid line) and between t = 5.9to 
and t = 6.6to (during the development of a large scale m = 2 dis- 
turbance, dashed line). 

|l999^ . We have indeed shown that the initial transient steep- 
ens the surface density profile in such a way that it matches 
the analytical estimates. The Mdisc = lAf* is even more 
complex: the disc never really settles down and it is self- 
regulated only in a time-averaged sense. The amplitude of 
the spiral disturbances changes strongly with time, and re- 
current, short-lived, low-m spiral patterns develop in the 
disc (see Fig. 110^ . A similar behavio ur has been also ob- 
served in t he cooling simulation s by ISellwood fc Carlberd 
iwM ) and lLauehlin et alJ lll998h . 

Another interesting results of our simulations is that, 
even at the high mass ratios considered here, they ap- 
pear to be still consistent with the fragmentation criterion 
tcooi < 3n~^ . The above criterion had been previously tested 
either in local (Gammie 2001|), or in global, low disc mass 
simulations ([I^i^e et al.n2 003al . Even if the global simula- 
tions bv iRlce et al .1 (|2003b^ show that the actual threshold 
value for the cooling time is slightly increased with respect to 
the local estimates by Gammic (2001), we have shown here 
that up to Mdisc = IMt, the condition tcooi > 7.557~^ effec- 
tively prevents the fragmentation of the disc. In fact, note 
that, even if rather large density perturbations are experi- 
enced by the disc, especially during the initial transients, 
none of these density enhancements is long-lived. This is 
due primarily to our choice for the cooling time scale and 
not to a limited resolution in our code. Indeed, simulations 
performed with the same number of particles as ours, but 
with a shorter cooling time, did show effective fragmentation 
(Rice ot al. 2003b). 

The present work completes the analysis described in 
Paper I, but can still be subject to many further refine- 
ments. In particular, all the results presented here and in 
Paper I are based on a specific form of the cooling term in 
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the energy equation. On the other hand, as aheady men- 
tioned, the outcome of gravitational instabilities in discs 
strongly depends on the assumptions about the disc ther- 
mal behaviour. Simulations that employ a cooling time con- 
stant with radius seem to show significant global energy 
transport iMeiia et alJl200^ . in contrast to our results. Of 
course, both a constant cooling timescale and our choice 
are highly idealized assumptions, and further studies are 
needed to assess which of the two is more reasonabl e. In 
principle, one could follow Ijohnson fc Gammid i2003h and 
find a relation between cooling time and density, based on 
some opacity prescription. This would not, in general, give 
tcooifi =const. However, a cooling time proportional to the 
dynamical timescale has the attractive feature of leading 
to a constant a, and since self-gravitating discs naturally 
evolve towards a self-regulated state, this might not be an 
unreasonable assumption (see also Bertin 1997}. 

Another important issue to be explored is how the disc 
response changes when different values for the ratio of the 
specific heats 7 is used and whether this can affect the frag- 
mentation criterion. Preliminary tests (Rice & Lodato, in 
preparation) seem to indicate that, when 7 is decreased, the 
boundary cooling time for fragmentation increases. Actually, 
it might be more reasonable to express the fragmentation 
criterion in terms of a maximum gravitational stress sus- 
tainable by a non-fragmenting disc. Fig. 1141 shows the rela- 
tionship between a and tcooi implied by Eq. Q , for different 
values of 7. The dotted line indicates the value of a corre- 
sponding to tcooi = 3n~^ when 7 = 2 (the value adopted by 
iGammielliiol . In Paper I and here, we have shown that, 
when the disc does not fragment, the strength of the gravi- 
tational disturbances is such that Eq. ||5J is approximately 
satisfied. If indeed there is a maximum a sustainable by a 
self-gravitating disc, then the disc would fragment whenever 
Eq. Q implies a larger value of a. This would then be con- 
sistent with an increase of the limiting cooling time when 7 
is decreased (see Fig. 1141 . Further investigations are clearly 
needed to asses whether this is the case or not. 
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